Giotto Suite Structure and Workflow
Jiaji George Chen
2023-08-02
1 Giotto Suite Object - A Flexible Spatial Analysis Framework
The Giotto object holds spatial-omic data and facilitates its analysis and visualization. Central to the design of the Giotto Suite object is the concept of spatial units and feature types. Most spatial data used to be pre-aggregated, with the expression information ‘belonging’ to an immutable unit of study. This lent themselves to similar analysis and organization as single-cell-based analyses. With recent methods, raw feature data is now decoupled from the morphological information, opening many doors for how the data can be analyzed and combined.
A spatial unit (spat_unit) is a set of polygonal information that can be used to define what you want your functional unit to be when performing analysis. This can be at different spatial sizes such as polygons that only select the nuclear region, the entire cell, or even larger tissue structures.
A feature type (feat_type) is which modality of feature you are analyzing such as rna or protein.
Analyses should then be organized first by which spat_unit is being studied, and then by which feat_type is being used to understand that spat_unit. Giotto’s updated object adopts a nested organization that reflects this in order to maximize flexiblity and ease of exploration.
2 Setup
i_p = installed.packages()
if(!"Giotto" %in% i_p)
devtools::install_github("drieslab/Giotto@suite")
# Ensure Giotto Data is installed
if(!"GiottoData" %in% i_p)
devtools::install_github("drieslab/GiottoData")
library(data.table)
library(Giotto)
library(GiottoData)
# Ensure the Python environment for Giotto has been installed
genv_exists = checkGiottoEnvironment()
if(!genv_exists){
# The following command need only be run once to install the Giotto environment
installGiottoEnvironment()
}3 Creating a Giotto Object
Assembling a complete object for spatial analysis requires expression
and spatial information.
This can be provided in the form of aggregated information (expression
matrices with paired spatial locations) or as raw feature data (feature
detections or images) and polygonal information.
This tutorial describes how you can assemble a Giotto object starting from subcellular data using a subset of the Vizgen MERSCOPE Mouse Brain Receptor Map data release.
3.1 Load Data
## provide path to vizgen folder
data_path = system.file('/Mini_datasets/Vizgen/Raw/',
package = 'GiottoData')
## 0.1 path to transcripts ####
# --------------------------- #
## each transcript has x, y and z coordinate
tx_path = paste0(data_path, '/', 'vizgen_transcripts.gz')
tx_dt = data.table::fread(tx_path)
tx_dt = tx_dt[,.(global_x, -global_y, gene, global_z)]
## 0.2 path to cell boundaries folder ####
# -------------------------------------- #
## vizgen already provides segmentation information in .hdf5 files
## Here I have already converted the hdf5 files to a simple data.table format
boundary_path = paste0(data_path, '/cell_boundaries/')
z0_polygon_DT = fread(paste0(boundary_path, '/', 'z0_polygons.gz'))
z1_polygon_DT = fread(paste0(boundary_path, '/', 'z1_polygons.gz'))
## 0.3 path to images ####
# ---------------------- #
DAPI_z0_image_path = paste0(data_path, '/', 'images/mini_dataset_dapi_z0.jpg')
DAPI_z1_image_path = paste0(data_path, '/', 'images/mini_dataset_dapi_z1.jpg')3.2 New Object Types
3.2.1 giottoPolygon
# giottoPolygon creation from data.table
z0_polygons = createGiottoPolygonsFromDfr(name = 'z0',
segmdfr = z0_polygon_DT)
z1_polygons = createGiottoPolygonsFromDfr(name = 'z1',
segmdfr = z1_polygon_DT)
z0_polygons## An object of class giottoPolygon
## spat_unit : "z0"
## Spatial Information:
## class : SpatVector
## geometry : polygons
## dimensions : 498, 1 (geometries, attributes)
## extent : 6399.244, 6903.243, -5152.39, -4694.868 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : poly_ID
## type : <chr>
## values : 40951783403982682273285375368232495429
## 240649020551054330404932383065726870513
## 274176126496863898679934791272921588227
## centroids : NULL
## overlaps : NULLplot(z0_polygons)
3.2.2 giottoPoints
# giottoPoints creation from data.table
tx = createGiottoPoints(tx_dt, feat_type = 'rna')
tx## An object of class giottoPoints
## feat_type : "rna"
## Feature Information:
## class : SpatVector
## geometry : points
## dimensions : 80343, 3 (geometries, attributes)
## extent : 6400.037, 6900.032, -5149.983, -4699.979 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : feat_ID global_z feat_ID_uniq
## type : <chr> <int> <int>
## values : Mlc1 0 1
## Gprc5b 0 2
## Gfap 0 3plot(tx)
3.2.3 giottoLargeImage
ultra_mini_extent = ext(c(6400.029, 6900.037, -5150.007, -4699.967 ))
image_paths = c(DAPI_z0_image_path, DAPI_z1_image_path)
image_names = c('dapi_z0', 'dapi_z1')
imagelist = createGiottoLargeImageList(raster_objects = image_paths,
names = image_names,
extent = ultra_mini_extent)
imagelist[[1]]## An object of class giottoLargeImage : "dapi_z0"
## Image extent : 6400.029, 6900.037, -5150.007, -4699.967 (xmin, xmax, ymin, ymax)
## Original image extent : 6400.029, 6900.037, -5150.007, -4699.967 (xmin, xmax, ymin, ymax)
## Scale factor : 0.108626547903541, 0.10862659908279 (x, y)
## Resolution : 9.2058527063567, 9.20584836903386 (x, y)
## Layers : 1
## Estimated max intensity : 208
## Estimated min intensity : 0
## Values : integers
## File path : '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/GiottoData//Mini_datasets/Vizgen/Raw//images/mini_dataset_dapi_z0.jpg'plot(imagelist[[1]])
3.3 Object Creation
Giotto objects are created through a call to
createGiottoObject(). Inputs to to this function should be
supplied as named lists of data where the names are used to determine
the data nesting they should be assigned as.
viz = createGiottoObject(feat_info = list('rna' = tx_dt),
spatial_info = list('z0' = z0_polygons))
instructions(viz, 'return_plot') = FALSE # set instructions to prevent returning of plots4 Subcellular Workflow
The next steps walk through the usual method of processing a Giotto object that has been assembled using subcellular information
4.1 Add Spatial Locations
Generate centroids to use as convenient spatial representations of
where your cells are. This information gets stored in the
@spatial_locs slot
viz = addSpatialCentroidLocations(gobject = viz,
poly_info = 'z0',
return_gobject = TRUE)
# visualize
spatPlot2D(viz)
4.2 Add Images
viz = addGiottoImage(viz, largeImages = imagelist[1])## Warning in spatInSituPlotPoints(viz, show_image = TRUE, largeImage_name = "dapi_z0", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)## select image done## plot image layer done## plot polygon layer done
4.3 Aggregate Expression Matrix
Using the polygonal information, sum up all of the transcript point
detections that fall within. This is can then be used as an expression
count matrix. This information is then stored in the
@expression slot.
# Calculate transcripts in 'rna' overlapped
# by the specified polygons in 'z0' (aggregation)
viz = calculateOverlapRaster(viz,
spatial_info = 'z0',
feat_info = 'rna')
# Send overlaps information to the expression slot as a new expression matrix
# of spat_unit 'z0' and feat_type 'rna' with a name of 'raw'
viz = overlapToMatrix(viz,
poly_info = 'z0',
feat_info = 'rna',
name = 'raw')## Warning in spatInSituPlotPoints(viz, show_polygon = TRUE, polygon_feat_type = "z0", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)## plot polygon layer done
From here on, the normal aggregate type analyses can continue.
4.4 Exploring the Object
We have added several objects to the Giotto object now. Directly
returning the object gives a summary of what information is contained
within.
For a more detailed look at a specific slot and type of data, Giotto has
several show functions.
# Returning the object provides a summary
viz## An object of class giotto
## >Active spat_unit: z0
## >Active feat_type: rna
## [SUBCELLULAR INFO]
## polygons : z0
## features : rna
## [AGGREGATE INFO]
## expression -----------------------
## [z0][rna] raw
## spatial locations ----------------
## [z0] raw
## attached images ------------------
## giottoLargeImage : dapi_z0
##
##
## Use objHistory() to see steps and params used# Closer look at @expression slot contents
showGiottoExpression(viz)## └──Spatial unit "z0"
## └──Feature type "rna"
## └──Expression data "raw" values:
## An object of class exprObj : "raw"
## spat_unit : "z0"
## feat_type : "rna"
## provenance: z0
##
## contains:
## 559 x 498 sparse Matrix of class "dgCMatrix"
##
## Mlc1 . . . . . . . . 1 . . 1 . ......
## Gprc5b . . 1 . 3 . . . 2 . 4 2 . ......
## Gfap . . . 2 1 1 . . 48 . 1 1 . ......
##
## ........suppressing 485 columns and 553 rows
##
## Chat . . . . . . . . . . . . . ......
## Blank-58 . . . . . . . . . . . . . ......
## Xcr1 . . . . . . . . . . . . . ......
##
## First four colnames:
## 40951783403982682273285375368232495429
## 240649020551054330404932383065726870513
## 274176126496863898679934791272921588227
## 323754550002953984063006506310071917306
## 5 Multiple Z-Stacks and Spatial Unit Provenance
What about when the
data is presented as multiple confocal z-stacks? The Vizgen MERSCOPE
Mouse Brain Receptor Map dataset is actually one such example, with 7
layers of different polygons alongside matching transcript detections at
each Z layer. This data should not be treated as fully 3D since the
z-layers are close enough together that they sample multiple times from
each individual cell as opposed to entirely different cells in a
different z layer.
This problem illustrates another way in which spat_unit
can be used in Giotto. Giotto handles this by defining each z-layer as
its own spat_unit and then performing separate
aggregations at each z-layer. These spat_units are then
merged them into a single combined spat_unit for
analysis.
Since the resulting combined spat_unit is different
depending on which layers were used, that information is tracked as
provenance.
5.1 Remake Giotto Object
In light of this, let’s remake our object.
# create with 2 spatial units
viz = createGiottoObject(feat_info = list('rna' = tx),
spatial_info = list('z0' = z0_polygons,
'z1' = z1_polygons))
instructions(viz, 'return_plot') = FALSE # set instructions to prevent returning of plots
# generate spatial locations
viz = addSpatialCentroidLocations(gobject = viz,
poly_info = paste0('z',0:1),
provenance = list('z0', 'z1'),
return_gobject = TRUE)5.2 Redo Aggregation
5.2.1 Layer-Specific Aggregations
viz = calculateOverlapRaster(viz,
spatial_info = 'z0',
feat_info = 'rna',
feat_subset_column = 'global_z',
feat_subset_ids = 0)
viz = overlapToMatrix(viz,
poly_info = 'z0',
feat_info = 'rna',
name = 'raw')
viz = calculateOverlapRaster(viz,
spatial_info = 'z1',
feat_info = 'rna',
feat_subset_column = 'global_z',
feat_subset_ids = 1)
viz = overlapToMatrix(viz,
poly_info = 'z1',
feat_info = 'rna',
name = 'raw')5.2.2 Collapse Into Common ‘aggregate’ spat_unit
# combine the calculated data for z1 and z0 into a new spatial unit called
# 'aggregate'
viz = aggregateStacks(gobject = viz,
spat_units = c('z0', 'z1'),
feat_type = 'rna',
values = 'raw',
summarize_expression = 'sum',
summarize_locations = 'mean',
new_spat_unit = 'aggregate')## Warning in spatInSituPlotPoints(viz, show_polygon = TRUE, polygon_feat_type = "aggregate", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)## plot polygon layer done
From here the usual aggregate analyses can be performed after filtering and normalization
viz = filterGiotto(viz,
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)
viz = normalizeGiotto(viz)##
## first scale feats and then cells# You can see that there are now z0, z1, and 'aggregate' as sets of spatial units
viz## An object of class giotto
## >Active spat_unit: aggregate
## >Active feat_type: rna
## [SUBCELLULAR INFO]
## polygons : z0 z1 aggregate
## features : rna
## [AGGREGATE INFO]
## expression -----------------------
## [z0][rna] raw
## [z1][rna] raw
## [aggregate][rna] raw normalized scaled
## spatial locations ----------------
## [z0] raw
## [z1] raw
## [aggregate] raw
##
##
## Use objHistory() to see steps and params usedshowGiottoExpression(viz)## ├──Spatial unit "z0"
## │ └──Feature type "rna"
## │ └──Expression data "raw" values:
## │ An object of class exprObj : "raw"
## │ spat_unit : "z0"
## │ feat_type : "rna"
## │ provenance: z0
## │
## │ contains:
## │ 494 x 498 sparse Matrix of class "dgCMatrix"
## │
## │ Mlc1 . . . . . . . . . . . . . ......
## │ Gprc5b . . 1 . 1 . . . 1 . 2 . . ......
## │ Gfap . . . 1 1 . . . 27 . . . . ......
## │
## │ ........suppressing 485 columns and 488 rows
## │
## │ Qrfpr . . . . . . . . . . . . . ......
## │ Chat . . . . . . . . . . . . . ......
## │ Blank-58 . . . . . . . . . . . . . ......
## │
## │ First four colnames:
## │ 40951783403982682273285375368232495429
## │ 240649020551054330404932383065726870513
## │ 274176126496863898679934791272921588227
## │ 323754550002953984063006506310071917306
## │
## ├──Spatial unit "z1"
## │ └──Feature type "rna"
## │ └──Expression data "raw" values:
## │ An object of class exprObj : "raw"
## │ spat_unit : "z1"
## │ feat_type : "rna"
## │ provenance: z1
## │
## │ contains:
## │ 494 x 504 sparse Matrix of class "dgCMatrix"
## │
## │ Mlc1 . . . . . . . . . . . 1 . ......
## │ Gprc5b . . . . . . . 2 . . . 1 . ......
## │ Gfap . 3 . 2 . 1 4 . . . . 18 . ......
## │
## │ ........suppressing 491 columns and 488 rows
## │
## │ Qrfpr . . . . . . . . . . . . . ......
## │ Chat . . . . . . . . . . . . . ......
## │ Blank-58 . . . . . . . . . . . . . ......
## │
## │ First four colnames:
## │ 40951783403982682273285375368232495429
## │ 17685062374745280598492217386845129350
## │ 223553142498364321238189328942498473503
## │ 240649020551054330404932383065726870513
## │
## └──Spatial unit "aggregate"
## └──Feature type "rna"
## ├──Expression data "raw" values:
## │ An object of class exprObj : "raw"
## │ spat_unit : "aggregate"
## │ feat_type : "rna"
## │ provenance: z0 z1
## │
## │ contains:
## │ 494 x 478 sparse Matrix of class "dgCMatrix"
## │
## │ Mlc1 . . . . . . 1 . . 1 . . . ......
## │ Gprc5b . 1 . 3 . . 2 . 4 2 . 1 . ......
## │ Gfap 2 . 2 1 . . 45 . 1 1 . . . ......
## │
## │ ........suppressing 465 columns and 488 rows
## │
## │ Qrfpr . . . . . . . . . . . . . ......
## │ Chat . . . . . . . . . . . . . ......
## │ Blank-58 . . . . . . . . . . . . . ......
## │
## │ First four colnames:
## │ 240649020551054330404932383065726870513
## │ 274176126496863898679934791272921588227
## │ 323754550002953984063006506310071917306
## │ 87260224659312905497866017323180367450
## │
## ├──Expression data "normalized" values:
## │ An object of class exprObj : "normalized"
## │ spat_unit : "aggregate"
## │ feat_type : "rna"
## │ provenance: z0 z1
## │
## │ contains:
## │ 494 x 478 sparse Matrix of class "dgCMatrix"
## │
## │ Mlc1 . . . . . . 5.537748 . . 6.080373 . . .
## │ Gprc5b . 6.658211 . 7.699722 . . 6.522136 . 7.927329 7.069674 . 5.338912 .
## │ Gfap 10.74423 . 8.16347 6.128572 . . 10.998911 . 5.945000 6.080373 . . .
## │
## │ Mlc1 ......
## │ Gprc5b ......
## │ Gfap ......
## │
## │ ........suppressing 465 columns and 488 rows
## │
## │ Qrfpr . . . . . . . . . . . . . ......
## │ Chat . . . . . . . . . . . . . ......
## │ Blank-58 . . . . . . . . . . . . . ......
## │
## │ First four colnames:
## │ 240649020551054330404932383065726870513
## │ 274176126496863898679934791272921588227
## │ 323754550002953984063006506310071917306
## │ 87260224659312905497866017323180367450
## │
## └──Expression data "scaled" values:
## An object of class exprObj
## for spatial unit: "aggregate" and feature type: "rna"
## Provenance: z0 z1
##
## contains:
## 494 x 478 dense matrix of class "dgeMatrix"
##
## [,1] [,2] [,3] [,4]
## Mlc1 -0.7655816 -0.5601177 -0.4625350 -0.4286685
## Gprc5b -1.3226931 2.0035912 -0.6914797 1.2966556
## Gfap 6.2252668 -0.8796777 1.6556304 0.6087384
## Ednrb 10.9380415 4.5135260 3.0888286 1.6443596
##
## First four colnames:
## 240649020551054330404932383065726870513
## 274176126496863898679934791272921588227
## 323754550002953984063006506310071917306
## 87260224659312905497866017323180367450
## 6 Autocorrelation
Giotto Suite also now includes geospatial statistics and can perform global and local spatial autocorrelation
6.1 Global Statistics
viz = createSpatialNetwork(viz,
method = 'kNN',
name = 'kNN_network',
k = 8)
# returns as a data.table of results
res = spatialAutoCorGlobal(
viz,
method = 'moran',
spatial_network_to_use = 'kNN_network',
data_to_use = 'expression',
expression_values = 'normalized'
)## No spatial weight matrix found in selected spatial network
## Generating distance matrix from kNN_networkprint((res))## feat_ID moran
## 1: Mlc1 6.278703e-02
## 2: Gprc5b 8.708937e-02
## 3: Gfap 3.521745e-01
## 4: Ednrb 1.488701e-01
## 5: Sox9 1.229596e-01
## ---
## 490: Blank-12 -2.874052e-03
## 491: Ffar4 -1.280580e-03
## 492: Qrfpr -3.748879e-03
## 493: Chat -3.093600e-03
## 494: Blank-58 -7.066559e-05print(res[moran == max(moran)])## feat_ID moran
## 1: Slc17a7 0.58388036.2 Local Statistics
# local autocorrelation results are returned as a set of spatial enrichment data and saved in the giotto object
viz = spatialAutoCorLocal(viz,
method = 'moran',
data_to_use = 'expression',
expression_values = 'normalized')## No spatial weight matrix found in selected spatial network
## Generating distance matrix from kNN_network## Attaching moran results as spatial enrichment: "expr_moran"# local moran's I of the most autocorrelated gene according to the global results
spatPlot2D(viz,
spat_enr_names = 'expr_moran',
cell_color = res[moran == max(moran), feat_ID],
color_as_factor = FALSE)
7 Generate Polygons
tx = getFeatureInfo(viz, feat_type = 'rna', return_giottoPoints = TRUE)
e = ext(tx) # get extent (spatial bounding box)
hexbin = tessellate(e, shape = 'hexagon', radius = 10, name = 'hexbin')
viz = setPolygonInfo(viz, hexbin)
pseudo_vis = makePseudoVisium(e)
viz = setPolygonInfo(viz, pseudo_vis, name = 'pseudo_visium')
viz = spatQueryGiottoPolygons(viz,
filters = list(pseudo_visium = 'all',
hexbin = 'all'))## Warning in spatQueryGiottoPolygons(viz, filters = list(pseudo_visium = "all", :
## NAs introduced by coercion
## Warning in spatQueryGiottoPolygons(viz, filters = list(pseudo_visium = "all", :
## NAs introduced by coercionplot(hexbin)
plot(pseudo_vis)
plot(viz@spatial_info$query_poly)
8 GiottoData Package
GiottoData is a set of mini Giotto objects that are already put together and have some analyses run on them. These were created as spatial subsets of the full data making them more portable and easy to manipulate.
mini_viz = GiottoData::loadGiottoMini('vizgen')
mini_viz9 Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] GiottoData_0.2.4 Giotto_3.3.1 data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.28 tidyselect_1.2.0 terra_1.7-39
## [4] xfun_0.39 bslib_0.4.2 listenv_0.9.0
## [7] lattice_0.20-45 colorspace_2.1-0 vctrs_0.6.2
## [10] generics_0.1.3 htmltools_0.5.5 yaml_2.3.7
## [13] utf8_1.2.3 rlang_1.1.1 R.oo_1.25.0
## [16] jquerylib_0.1.4 pillar_1.9.0 glue_1.6.2
## [19] withr_2.5.0 R.utils_2.12.2 lifecycle_1.0.3
## [22] progressr_0.13.0 munsell_0.5.0 gtable_0.3.3
## [25] R.methodsS3_1.8.2 future_1.32.0 codetools_0.2-18
## [28] evaluate_0.21 labeling_0.4.2 knitr_1.42
## [31] fastmap_1.1.1 parallel_4.2.1 fansi_1.0.4
## [34] highr_0.10 Rcpp_1.0.11 backports_1.4.1
## [37] checkmate_2.2.0 scales_1.2.1 cachem_1.0.8
## [40] dbscan_1.1-11 jsonlite_1.8.4 parallelly_1.35.0
## [43] farver_2.1.1 ggplot2_3.4.2 png_0.1-8
## [46] digest_0.6.31 dplyr_1.1.2 grid_4.2.1
## [49] rprojroot_2.0.3 scattermore_0.8 here_1.0.1
## [52] cli_3.6.1 tools_4.2.1 magrittr_2.0.3
## [55] sass_0.4.6 tibble_3.2.1 future.apply_1.10.0
## [58] pkgconfig_2.0.3 Matrix_1.5-4 rmarkdown_2.21
## [61] rstudioapi_0.14 globals_0.16.2 R6_2.5.1
## [64] igraph_1.4.2 compiler_4.2.1